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We present a model based on the lattice Boltzmann equation that is suitable for the simulation 
of dynamic wetting. The model is capable of exhibiting fundamental interfacial phenomena such as 
weak adsorption of fluid on the solid substrate and the presence of a thin surface film within which 
a disjoining pressure acts. Dynamics in this surface film, tightly coupled with hydrodynamics in 
the fluid bulk, determine macroscopic properties of primary interest: the hydrodynamic slip; the 
equilibrium contact angle; and the static and dynamic hysteresis of the contact angles. The pseudo- 
potentials employed for fluid-solid interactions are composed of a repulsive core and an attractive 
tail that can be independently adjusted. This enables effective modification of the functional form 
of the disjoining pressure so that one can vary the static and dynamic hysteresis on surfaces that 
exhibit the same equilibrium contact angle. The modeled solid-fluid interface is diffuse, represented 
by a wall probability function which ultimately controls the momentum exchange between solid and 
fluid phases. This approach allows us to effectively vary the slip length for a given wettability (i.e. 
the static contact angle) of the solid substrate. 



PACS numbers: 47.85.-g;47.55.-t;68.08.Bc 

I. INTRODUCTION 

Dynamic wetting encompasses a class of interfacial 
phenomena relevant to technical applications in strategic 
areas ranging from materials science to energy research. 
Critical applications such as coating, self-assembly, and 
microfluidic handling have motivated extensive efforts to- 
ward developing models for the dynamic wetting of solid 
surfaces. Nonetheless, physical understanding and suc- 
cessful mathematical modeling of these phenomena still 
requires significant developments l2|] . A comprehensive 
description of dynamic wetting must account for the cou- 
pling between molecular and hydrodynamic interactions 
taking place on length and time scales that extend over 
several orders of magnitude. 

Classical continuum-level descriptions rely on the 
Navier- Stokes equations for the mathematical modeling 
of hydrodynamics, while microscopic interactions need 
to be coarse grained in order to render hydrodynamic 
boundary conditions at fluid-solid or fluid-fluid inter- 
faces. In order to model the physical coupHng between 
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microscopic interactions and hydrodynamics, boundary 
conditions for the wall velocity and the static/dynamic 
contact angles commonly take the form of algebraic equa- 
tions [3] or even partial differential equations [4] . From a 
mathematical viewpoint, these boundary conditions are 
not just input parameter values, but rather involve addi- 
tional model equations. Furthermore, the set of hydro- 
dynamic equations must usually be supplemented with 
equations of state (EOS) and other constitutive relations 
that can make a continuum-level description cumber- 
some. 

This type of classical continuum description, although 
widely employed, suffers from several limitations in deal- 
ing with nontrivial microscopic effects. In particular, the 
sharp-interface limit leads to well-known mathematical 
singularities and paradoxes; e.g. logarithmic stress sin- 
gularities, the moving contact line paradox [El, [6|. Such 
unphysical artifacts are readily removed [g] by consider- 
ing that fluid-fluid and fluid-solid interfaces are actually 
"diffuse" and thus have a finite thickness, determined 
by thermal diffusion and the range of action of molecu- 
lar forces. The interface thickness can be considered as 
a cutoff length [3| , or characteristic microscopic scale of 
the physical system, below which the classical continuum 
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model with a sharp interface is not vahd. 

A fundamental phenomenon observed at the micro- 
scopic scales is the presence of thin interfacial films (e.g. 
precursor films on hydrophilic surfaces [6|, |7|) within 
which dispersion forces (e.g. London-van der Waals 
forces) become significant. The mass and momentum 
fluxes through these interfacial films can produce non- 
trivial effects on the hydrodynamic behavior (e.g. effec- 
tive slip velocity 0,0,111) and the equilibrium conditions 
prevailing in the fluid bulk [q] . It becomes therefore nec- 
essary to refine the level of modeling within the interfacial 
films in order to model dynamic wetting phenomena in 
many systems of practical interest (e.g. micro /nanoflows, 
colloids). 

In this work, we propose a so-called "mesoscopic" 
model fld-^ embodying an augmented description of 
the soHd-fluid interface. Such a mesoscopic approach, 
based on the single-particle distribution, is convenient 
for modeling the net effects of microscopic interaction 
forces on the dynamics of macroscopic quantities (e.g. 
mass density, thermodynamic pressure, and fluid mo- 
menta). This strategy can be realized by using mean- 
field interaction potentials, or pseudo-potentials, that are 
scalar functions of macroscopic quantities (e.g. mass 
density). Similarly to effective interaction potentials 
used in different forms of DLVO theory [7, 9], pseudo- 
potentials are scaled by interaction parameters (e.g. at- 
traction/repulsion Hamaker constants). These interac- 
tion parameters determine an interfacial force per unit 
area, or disjoining pressure [H, [l3|, [l^ , that can have both 
repulsive and attractive components varying as a func- 
tion of the distance from the interface. The disjoining 
pressure near the solid-fluid interface ultimately deter- 
mines the contact angle at the apparent three-phase con- 
tact line [3,|9|, Ell- Hence, the set of attraction/repulsion 
parameters employed determines implicitly the wetting 
properties of the solid. Moreover, the attractive-repulsive 
character of the modeled interactions can lead to multi- 
ple local minima and maxima in the disjoining pressure. 
These physical features give rise to nontrivial effects, such 
as the static hysteresis of the apparent contact angle on 
a (macroscopically) smooth surface [9]. 

This article is structured as follows. In Sec. II we 
formulate our mesoscopic model motivated by the physi- 
cal insights presented above; the current implementation 
of this model is based on the numerical solution of a lat- 
tice Boltzmann (LB) equation. A fundamental difference 
with previous LB models lies in the augmented treatment 
of the solid-fluid interface by means of a wall probability 
function and a pseudo-potential for solid-fluid interac- 
tions having two adjustable parameters that control the 
magnitude of attractive and repulsive interfacial forces. 
The two-parameter model allows us to simulate surfaces 
that exhibit the same static properties, i.e. same static 
contact angles, yet different dynamic behavior, i.e. dif- 
ferent static or dynamic contact angle hystereses. In Sec. 
Ill we report numerical results for pressure-driven flows 
of a volatile fluid. We focus on the static contact angles 



of sessile droplets and advancing/receding contact angles 
of droplets under pressure-driven flow. The results reveal 
that the proposed model reproduces key features of real 
soHd surfaces. In Sec. IV we conclude with a summary 
of the key results and outline potential directions for the 
presented approach. 



II. A MESOSCOPIC MODEL FOR 
HYDRODYNAMICS AND INTERFACIAL 
PHENOMENA 

The model we propose is based on the Boltzmann- 
BGK equation for a single-component fluid [lo|, [HI 
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Here / = f{-K^v,t) is the single-particle probability dis- 
tribution at position coordinate x and velocity coordinate 
V at a time instance t. The term Sf/6t on the right-hand 
side of Eq. [1] accounts for the action of both external and 
internal forces. The single relaxation time r is a model 
parameter that determines the kinematic viscosity, the 
mobility, and other transport coefficients. The equilib- 
rium distribution function is a Maxwell-Boltzmann dis- 
tribution 
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where D is the velocity space dimensionality, p is the 
local mass density, and u is the fluid velocity. The dis- 
tribution variance = kBT/m is determined by the spe- 
cific thermal energy (T: fluid temperature, ks'- Boltz- 
mann constant, m: molecular mass) and defines the 
characteristic scale of diffusive processes. For the stud- 
ied single-component system we adopt a unit molecular 
mass, m = 1, without loss of generality. The LB model in 
this work is valid for isothermal flows (T=const.) where 
a classical hydrodynamic description involves the three 
leading moments of the distribution function: 



M(0)(x,t) = j /(x,v,t)dv 



p; 



M(i)(x,t) = j /(x,v,t)vdv = pu; 



M^2)(x,t) = y'/(x,v,t)vvdv : 



puu 
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The second-order moment, M^^\ is the momentum flux 
tensor given by the sum of the macroscopic momentum 
flux, puu, and the stress tensor, tr; in a closed-form de- 
scription, cr is a functional of mass density and fluid ve- 
locity. A comparable level of hydrodynamic description 
is attained with a second-order LB method [ll,il5,] where 
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Vz or Vi ^ states Wi 
(±1,0),(0,±1) 1-4 1/12 

(±1,±1) 5-8 2/27 

(±2,0), (±2,0) 9-12 7/360 

(±2, ±2) 13-16 1/432 

(±3,0),(0,±3) 17-20 1/1620 
(0,0) 21 91/324 

" = mv^ = 2/3 

TABLE I. Lattice D2Q21. Velocity abscissas (in lattice units) 
and Gauss-Hermite weights. 

the distribution function / is a truncated Hermite expan- 
sion 

/(x,v,t) = /^(v)[M(^)±iM(i):v (6) 
+ (M(2) - M^^^Ol) : (vv - 01)]. 

Here = {27t0)-^^^ exp(-vV2(9) is a Gaussian weight 
and I is the unit tensor. High-order moments (M^^^; n > 
2) are readily evaluated in terms of the three leading or- 
der moments (M^^^; n = 0, 2) by using Eq. [6l A projec- 
tion step 0,|T5| needs to be included in the LB algorithm 
in order to enforce the functional form in Eq. [6l 



A. Lattice Implementation 

The numerical algorithm evolves a finite set of distri- 
butions /i(x, t) = /(x, v^,t) {i = 1,(5). The lattice ve- 
locities (v^; i = lj>Q) are defined by a Gauss-Hermite 
quadrature rule [IJJ. The quadrature rule must fulfill a 
fourth-order algebraic degree of precision {d > 4) to ac- 
curately approximate the isothermal flow solution in the 
continuum limit. In the presence of large density gradi- 
ents that can develop at an interface, a high-order rota- 
tional symmetry is highly desirable after discretization 
of velocity space. Such symmetry is critical to effectively 
retain the isotropy of high-order spatial derivatives re- 
quired for the discrete approximation of interaction forces 
[l6j . For the current implementation we adopt a D2Q21 
lattice (dimensions D = 2 and states Q = 21)pjj; ab- 
scissas and corresponding weights are reported in Tab. HI 
The D2Q21 lattice velocities are the integration points 
of a quadrature rule having seventh-order algebraic pre- 
cision and satisfy moment isotropy up to the sixth-order 
[ll]. This lattice retains isotropy of the fifth-order spa- 
tial derivatives in the discrete gradient operator (Eqs. [7]- 
[8]) [16]. Thus, the D2Q21 lattice has comparatively better 
performance over low order lattices (e.g. D2Q9, D3Q17) 
in terms of reducing spurious currents and other numer- 
ical artifacts. The main disadvantage of using a high 
order lattice is that numerical simulation of even sim- 
ple geometries can become computationally intensive. A 
lattice cell (showed in Fig. [1]) extends over six nodes; 
i.e. each node propagates information up to three lattice 




FIG. 1. D2Q21 lattice cell. This lattice satisfies sixth-order 
moment isotropy. Discrete gradient operators implemented 
on D2Q21 lattices retain isotropy up to the fifth-order spatial 
derivatives. 



sites away in one time step. With the implementation in 
this work, properly resolving an interface or effectively 
rendering a boundary condition requires a minimum res- 
olution of one lattice cell (i.e. six nodes). For simulating 
nano- or microscale flows the computational cost of the 
present approach remains affordable by current compu- 
tational resources relatively low when compared 
to molecular dynamics simulations. Among other ad- 
vantages of the proposed approach, the LB algorithm is 
particularly suitable for massive parallelization and /or 
acceleration using graphic processing units (GPUs) jlSj . 



B. Non-ideal fluid behavior and interfacial 
phenomena 

The mesoscopic description based on Eq. [1] employs 
the BGK ansatz to model the macroscopic effect of short- 
range fluid- fluid interactions. All other interaction forces, 
responsible for non-ideal fluid behavior (e.g. non-ideal 
equation of state, phase separation) and interfacial phe- 
nomena (e.g. surface tension, disjoining pressure, par- 
tial wetting) are modeled by an approximation of the 
actual force term df/St = —m~^g - df /dv in the ki- 
netic transport equation; g is a body force that may 
also include external fields. Among different possible 
approximations for Sf / St [III, [HI we adopted the exact 
difference method introduced in [19^ because of its op- 
timal numerical stability [20]. Hence, the force term in 
Eq.[T]is Sf /St = /^^(p, u*) — /^^(p, u) where the equilib- 
rium distribution is computed using a "shifted" velocity 
u* = u + gAt. 

For the class of LB models we employ [ll|, |l2|, the 
volumetric body force pg = ?/^(x, t)V / i(;(|r|)?/;(x + 
r, t)dr is determined by a spatial convolution of the 
pseudo-potentials with the Gaussian kernel i^drl) = 
(27rA>:l9)-2/^exp(-|rp/2/^6>). The characteristic length 
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scale of the interaction kernel determines the thick- 
ness of the resulting interfaces (in all our computations 
we use = 1). The interaction force pg = Fff -\-Ffs 
contains a Fluid-Fluid component 

Q 

Fff (x, t) = ipFF (x, t) ^ w^ipFF (x + , t)Ti , (7) 

and a Fluid-Solid contribution 

Q 

FFsi^^t) = p{x,t)'^WilljFs{^^T^i^t)Ti^ AFsi^.t). 

i=0 

(8) 

The long-range fluid-solid interactions thus receive sim- 
ilar treatment to that described in [l2|, using different 
pseudo-potentials for cross interactions between species. 
Hereafter, we refer to all those microscopic interactions 
that act beyond an atomic radius as long-range; in this 
work long-range forces are mainly associated to van der 
Waals interactions and ionic double layers. 

The last term, AF^, introduces a momentum ex- 
change between fluid and solid molecules attributed to 
short-range interactions (considered as probabilistic col- 
lision events) occurring in the region adjacent to the solid 
surface. The discrete lattice directions, r^, and Gauss- 
Hermite weights, Wi^ are the same as those employed for 
the D2Q21 lattice [see Tab.l]. All dimensional quantities 
in this work are reported in lattice units. 

Fluid-Fluid interactions. The Fluid-Fluid interac- 
tions are defined as: 

^FF (x, t) = y/2[p0-pEos{p.0)i (9) 

where the pressure Peos is given by an equation of state 
(EOS). In this work, we model a volatile fluid that can 
separate into a (stable) vapor and liquid phase at the 
studied temperature T = mO/kB, while both phases ex- 
hibit ideal fluid behavior (i.e. p (x p). The adopted EOS 
is given by a piecewise linear relation [20[ 

{pOv P< Pi 

Pi^{p-pi)Ou Pi<P<p2 (10) 
P2 + (P - P2)0l P> P2 

where pi and p2 are the endpoints of the unstable branch. 
Pi = PiOv and p2 = piOy + (p2 - pi)Ou, and Oy > 0, 
Ou < 0, and 6l > are the slopes in the vapor, unstable, 
and liquid branches, respectively. The parameters em- 
ployed in the present work are: Oy = 0.25^; Ou = —0.25^; 
Ol = 1.0(9; pv = 0.1; pi = 0.222; p2 = 0.869; and 
PL = 1.0. Based on the continuum calculation of thermo- 
dynamic equilibrium, this parameter combination pro- 
duces phase equilibrium at a density ratio Pl/pv = 10 
and a compressibility ratio /3 = Pl^l/ Pv^v = 1/40. For 
this system we report a surface tension 7 c:^ 0.09 com- 
puted in simulations via force integration across a planar 
interface, as well as by measuring pressure differences 
across a circular interface [for a description of the proce- 
dure see 19, and • 



Fluid-Solid interactions. The solid phase is treated 
on a similar footing with the liquid phase. The proper 
choice of ipFS can effectively model a disjoining pres- 
sure acting within a surface film of finite thickness. The 
pseudo-potential employed for Fluid-Solid interactions 
has the general form 

^FsM = Gni^ni^) + Ga^a(x). (11) 

A procedure to determine the repulsive, i/jr^ and attrac- 
tive potentials, i/ja^ for arbitrary surface geometries is 
described in the Appendix. Once the pseudo-potentials 
are properly defined, the static contact angle Oy is mod- 
ified by adjusting the attraction parameters Gr and 
Ga' In analogy with DLVO theory, the modeled Fluid- 
Solid interactions have a repulsive component Gr^i/jr (e.g. 
attributed to double-layer repulsion) and an attractive 
component Ga^^a (e.g. due to London-van der Waals 
forces). The fundamental feature modeled in Eq. [TTI is 
that the resulting surface forces reverse direction as the 
solid is approached. The term in Eq.[8]that models short- 
range interactions (i.e. elastic collisions, Pauli repulsion), 

AFs(x,t) = p(x,t)(/)^(x)Auw^(x,t), (12) 

produces a velocity shift Au^aii = ^waii — u — Fff At 
in the fluid adjacent to the solid. The total momen- 
tum exchange between fluid and solid dynamically deter- 
mines the hydrodynamic slip. In Eq. [T2]we consider that 
the solid surface can present a microscopic scale rough- 
ness that results in a diffuse fluid-solid interface for the 
present mesoscopic description. Hence, the wall prob- 
ability function 05 (x) takes finite values < <ps < ^ 
within the interfacial region where thin surface films de- 
velop. While (j)s = I inside the solid bulk, the collision 
probability vanishes (^5- = 0) at a certain distance from 
the solid bulk. In practice, the spatial variation of (j)^ 
allows us to model the influence of microscale roughness 
on the effective slip [1]. Once the wall function ^5 is 
defined (e.g. using the procedure in the Appendix), the 
exponent e > 1 in Eq. [12] provides a means of incorpo- 
rating different levels of microscale roughness via "fine" 
adjustments to the interface sharpness. 

III. NUMERICAL RESULTS 

In this section, we present numerical results that 
demonstrate the capabilities of the method formulated 
in Sec. II and we discuss briefly the most relevant obser- 
vations. All simulations are performed employing small 
values of the relaxation time r = 1.0-1.5; within this 
range the reported results present no significant depen- 
dence on r. The only EOS employed and the correspond- 
ing equilibrium conditions are described in Sec. II. The 
functions 05, i/jr, and i[ja that determine the solid loca- 
tion and all Fluid- Solid interactions are evaluated at ini- 
tialization (i.e. before the actual dynamic simulation) us- 
ing the numerical procedure described in the Appendix. 
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FIG. 2. Wall functions determined before dynamic simula- 
tion: wall probability (e = 1, 1.3); repulsive and attractive 
components of V^fs = Gri/jr + GaiPa; y (in lattice units) is in 
the direction normal to a wall (bottom or top). The lines in 
the plots are interpolation curves fitting the computed node 
values [see Appendix]. The value of the wall probability is 
used to determine three regions: (I) solid substrate (^s > 0.5) 
where short-range (fluid-solid) interactions are dominant; (II) 
interfacial boundary layer (0.5 < < 0.01) where surface 
forces are significant; (III) hydrodynamic region or fluid bulk 
(05 < 0.01) where surface forces become neghgible. 



FIG. 3. Fluid- solid potential (Top panel), iprs, and volu- 
metric force (Bottom panel), Qw, for interaction parameters 
within the range employed in simulations: Gr = 1.33; Ga = 
0.2, 0.5, and 0.9. The lines in the plots are interpolation 
curves fitting the computed node values [see Appendix] . The 
^-axis is in the direction normal to a wall (y coordinates are 
indicated in lattice units). The hydrodynamic region (I) be- 
gins at yw — 8.5 where (psi^y-w) = 0.01 (the opposite boundary 
is located at y^ — Ly — 7.5). The normal component Qw (given 
in lattice units) gives rise to a disjoining pressure 11 = —pQw 
that dominates within the interfacial region (II). 



The wall functions employed for the present simulations 
are reported in Fig. [2l We present in Fig. [3] the Fluid- 
Solid potential iIjfs and the volumetric force g^^ for rep- 
resentative values of the interaction parameters. For the 
simulated flat surfaces only the normal component of g^^ 
is active, causing a disjoining pressure H = —pQw The 
curves in Fig. [3] illustrate the effects of scaling the attrac- 
tive interaction tail of ^|)FS] this ultimately controls the 
surface wettability. 

The modeled microscopic interactions give rise to three 
distinct spatial regions [see Figs. [2]-[3] where different 
effects dominate the dynamics. The value of the wall 
probability (0 < (/>s' < 1) determines the relative impor- 
tance of fluid-solid interactions near the solid, and their 
effective absence within the fluid bulk. We use charac- 
teristic values of the wall probability to define such re- 
gions: (I) solid substrate {(j)s > 0-5) where short-range 
interactions are dominant and a vapor fraction is ad- 
sorbed {p/pL ^ 0.001 for Or = 3.0, p/pL ^ 0.03 for 
Gr = 0.9); (II) interfacial film (0.5 < (j^s < 0.01), or 
boundary layer, where long-range interactions and the 
resulting disjoining pressure 11 = —pgw dominate; (III) 
fluid bulk {(j)s < 0.01), or hydrodynamic region, where 
surface forces are negligible and the Navier- Stokes equa- 
tions are recovered. In our simulations, the first node in- 
side the hydrodynamic region {(j)s ^ 0.01) is located at 9 
lattice units, or 1.5 lattice cells, from the outer boundary 
of the simulation domain. The hydrodynamic velocity on 
the wall is observed dit = 8.5, which corresponds to our 
limit {(ps — 0.01) for the interfacial region {y^ = Ly — 7.5 
for the opposite wall with Ly being the domain height). 



A. Pressure-Driven Flow in capillaries 

Before studying problems that involve contact lines, 
we simulate pressure-driven flow of a volatile liquid 
{pL = 1-0) in a two-dimensional channel of dimensions 
Lx X Ly = 10 X 74 with periodic boundary conditions in 
the X— direction. A small pressure difference Ap = dpxLx 
is applied in the x-direction via a small body force 
dpx = 4.0 X 10~^/p, while we vary the interaction pa- 
rameters that control surface wettability. In order to 
model the effects of varying the microscopic-scale rough- 
ness we adjust the exponent of the wall probability func- 
tion this exponent effectively adjusts the momentum 
exchange due to short-range repulsive interactions within 
the interfacial region (II) where long-range surface forces 
are active. 

The results in Figs. |4]-[5] report mass density and fluid 
momentum profiles for three sets of interaction parame- 
ters {Or = 1.33, Ga = 0.2, 0.5, 0.9). We compare two 
values of the wall function exponent: e = 1.0 to model a 
"rough" surface where microscopic roughness and long- 
range interactions have the same characteristic length; 
and e = 1.3 to model a "smooth" surface where the mi- 
croscopic roughness is smaller than the range of action of 
long-range surface forces. The density profile in Figs.|4]-[5] 
shows an interfacial film that develops between the hy- 
drodynamic region and the solid substrate. We observe 
that the thickness and density of the surface films are 
adjusted by varying the interaction parameters, Gr and 
Ga, that control the Fluid-Solid, tprs^ and Fluid-Fluid, 
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(a) Gr = 1.33, Ga = 0.2, e = 1.0 (a) Gr = 1.33, Ga = 0.2, e = 1.3 




(b) Gr = 1.33, Ga = 0.5, e = 1.0 (b) = 1.33, Ga = 0.5, e = 1.3 




(c) = 1.33, Ga = 0.9, e = 1.0 (c) = 1.33, Ga = 0.9, e = 1.3 



FIG. 4. Pressure- driven flow on a microscopically "rough" 
surface (e = 1.0) for different surface wettability conditions 
{Gr = 1.33, Ga = 0.2, 0.5, and 0.9). The (red) solid line 
indicates the analytical solution of incompressible Poiseuille 
flow, which is valid within the hydrodynamic region (III), 
for Ap = 4.0 X 10~^ and using a no-slip boundary condition 
u{yw) = 0. A surface film of varying density develops within 
the interfacial region (II). 



FIG. 5. Pressure-driven flow on a microscopically "smooth" 
surface (e = 1.3) for different surface wettability conditions 
{Gr = 1.33, Ga = 0.2, 0.5, and 0.9). The (red) solid line 
indicates the analytical solution of incompressible Poiseuille 
flow, which is valid within the hydrodynamic region (III), 
for Ap = 4.0 X 10~^ and using a no-slip boundary condition 
u{yw) - 0. 



tpFF^ pseudo-potentials. It follows that the modeled 
physico-chemical properties of the fluid and solid phases 
determine the hydrodynamic slip. The amount of hy- 
drodynamic slip is determined by comparing against the 
analytical solution of incompressible Poiseuille flow with 
no-slip velocity applied at i/uj . In agreement with molec- 



ular dynamic simulations of similar pressure-driven flows 
[8], the amount of effective slip is rather independent of 
the studied flow conditions. As seen in Figs.HHSl there is 
no clear correlation between the hydrodynamic slip and 
the surface wettability, determined by the attraction pa- 
rameter {Ga = 0.2-0.9) when the repulsion parameter 
is fixed {Gr = 1.33). The scale of the (modeled) micro- 
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roughness, controlled by e, has significant effects on the 
hydro dynamic slip. A microscopically "rough" surface 
(e = 1.0) presents no effective slip [Fig. |4], regardless of 
its wettability, while a "smooth" surface exhibits signifi- 
cant slip [Fig. [5] even for hydrophilic conditions produced 
by a high attraction parameter {Ga = 0.9). We remark 
that these qualitative physical features modeled with the 
present approach are reported in a large number of ex- 
perimental studies [1]. 

B. Static wetting and equilibrium contact angle 

We simulate a two-dimensional drop on a flat surface 
in order to quantify the surface wettability when vary- 
ing the surface interaction potentials through the pa- 
rameters Gr and Ga' The simulation domain size is 
Lx X Ly = 650 X 150 lattice units and the drop vol- 
ume per unit width is Vd = ttRq; the results present no 
significant dependence on the droplet volumes employed 
{Rq = 30-80) within the studied range of contact an- 
gles 20° < < 160°. The equilibrium contact angle 
is determined by numerically fitting a circle of radius R 
[see Figs. [6]-[7] to the vapor-liquid interface which is de- 
fined by the contour line for p = {pL Pv)/'^- The circle 
fitting by least squares is confined to the hydrodynamic 
region (^5 > 0.01), above the interfacial film, where sur- 
face forces are negligible and constant curvature of the 
droplet is to be expected [9]. The reported equilibrium 
contact angle Oy = acos(l — h/R) is evaluated using the 
circle radius R and droplet height h (i.e. distance be- 
tween the apex and the bottom liquid-vapor interface) 
[see Figs. [BHl]. The apparent contact angle, which we do 
not report in this work, should be measured above the 
interfacial surface film where y = y^. 

The density field reported in Figs. [6H3 exhibits funda- 
mental physical features observed in previous work ac- 
counting for the presence of van der Waals forces [see 
[H, [il, [13I, ]• The most relevant feature observed in 
our simulations is the development of interfacial films 
between the sessile droplet and the solid substrate; in 
Fig. [71 we identify the overlapping boundary regions an- 
alytically studied in |9l]. As expected, interfacial forces 
adjusted via Ga and Gr determine the position of the 
vapor-liquid interface {p = 0.55), at which mechanical 
equilibrium is attained, as well as the equilibrium contact 
angle. The equilibrium contact angle is reported in Fig. [8] 
as a function of the attraction parameter Ga for three dif- 
ferent values of the repulsive parameter Gr and the mod- 
eled micro-roughness effects studied in Sec lIII A l Emp loy- 
ing the augmented Young-Laplace equation [7|, Isl, [l^ the 
equilibrium contact angle is determined by 

1 

cos6>y = 1 + - / Ii{y)dy (13) 
7 Jo 

= -[Eo{GR)^kAGA]^0{G\) 

7 

after adopting a linear approximation for the surface en- 
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(d) 

FIG. 6. Sessile circular droplet at equilibrium. The equi- 
librium contact angle Oy — acos(l — h/R) is computed via 
circle fitting (solid black line); R is the circle radius, and h 
indicates the droplet height (measured from the outer bound- 
ary of the simulation domain). Contour lines (solid) for the 
fluid density, p = 0.01,0.1,0.55, indicate a gradual growth 
of the interfacial film thickness as the attraction parameter 
Ga increases. Contour lines (dashed) for the solid proba- 
bility function, p^ — 0.1,0.5, contain the interfacial region 
where a disjoining pressure is active. All cases correspond 
to Or = 1.33 and Rq = 50. (a) Ga = 0.0, By = 150°. 
(b) Ga = 0.2, Oy = 112°. (c) Ga = 0.5, Oy = 58.7°. (d) 
Ga = 0.7, Oy = 36.7°. 
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(a) 




FIG. 7. Sessile circular droplet on a hydrophilic surface. 
Circle fitting and contour line description is as in Fig. [6] 
for Gr = 1.33 and Rq = 50. (a) Top panel: Ga = 0.8, 
Oy = 30.9°; bottom panel: Ga = 0.9, By = 24.0°. (b) Den- 
sity field in the vicinity of the (apparent) three-phase con- 
tact line [see [9|]]. I: Solid bulk, where a small vapor fraction 
(p < 0.001) is adsorbed; II- 1: interfacial film at the liquid- 
solid interface; II- 2: overlap region between liquid- vapor and 
liquid-solid interfaces; II-3: thin interfacial film at the vapor- 
solid interface. Ill: liquid bulk. IV: liquid-vapor interfacial 
region. Interfacial forces (e.g. disjoining pressure) in regions 
II- (1-3) are determined by both interaction pseudo-potentials 
il^FS and V^FF- 



ergy difference "excess" = Eq{Gr) + A:aGa. The 
linear approximation in Eq. [13] fits the numerical data in 
Fig. [8fa-b) for moderate-to-large values of the contact 
angle. 

Similar behavior of the contact angle variation with 
the surface energy excess, and an apparent saturation of 
the minimum contact angle attained, has been reported 
by an alternative approach based on Cahn-Hilliard mod- 
els for the surface free-energy [H, [23|. The speculated 
reason for the apparent saturation of the minimum con- 
tact angle is geometric effects [23]; because the size of 
the droplet is comparatively small with respect to the 
surface film, at low but finite contact angles the droplet 
transitions into a thick film and the circle fitting pro- 
cedure becomes inaccurate. Nevertheless, the reported 
behavior at low contact angles presented no significant 
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FIG. 8. Equilibrium contact angle Oy for a two-dimensional 
droplet on a flat surface versus attraction parameter Ga ■ The 
reference radius Ro determines the droplet volume (per unit 
width) Vd = ttRI. The same equilibrium contact angle can 
be obtained using infinite combinations of wall interaction 
parameters Gr and Ga- Solid lines are estimations by an 
augmented Young- Laplace equation cosOy ~ {Eo-\-kAGA)/j' 
(I) ^0 = -0.64, kA = 2.5; (II) Eq = -0.87, /ca = 2.2; (III) 
Eo = -1.92, kA = 1.76; (IV) Eq = -1.0, kA = 3.06; in 
all cases 7 = 0.09. Deviations from a linear approximation 
for the surface energy excess AE = kAGA are observed for 
Oy < 45°. 



dependence on the droplet sizes employed {Rq = 30-80) 
in our simulations. Compressibility effects can also con- 
tribute to the observed deviations from a linear increase 
in the surface energy AE with respect to Ga and the ap- 
parent saturation to a minimum contact angle attained. 
A key observation in this section is that infinite combi- 
nations of interaction parameters (Gr^Ga) can produce 
the same equilibrium contact angle. In the following sub- 
section we study dynamic wetting properties of surfaces 
that exhibit the same static contact angle, attained by 
different combinations of interaction parameters. 
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h=55. 8 Advancing interface 
circle fit: R=109.4 




FIG. 9. Dynamic contact angles for a moving droplet. Ad- 
vancing/receding contact angles, Oa and Or^ are determined 
by fitting two circles, of radii i?, to the vapor-liquid inter- 
face from the apex, at height /i, to the corresponding ad- 
vancing/receding contact lines. The illustrated case corre- 
sponds to Ca = 0.006 using interaction parameters Gr = 1.33 
and Ga — 0.5, which produce an equilibrium contact angle 
By - 58°. 



C. Contact Angle Dynamics 

After characterizing static wetting properties of the 
modeled surfaces we investigate dynamic wetting con- 
ditions by applying a pressure difference that causes a 
flow in the x-direction and motion of the droplet; these 
simulations are performed for the same domain as in 
Sec. IIIIBl The same static contact angle Oy — 58° 
is obtained with three sets of interaction parameters: 
{Or = 1.33, Ga = 0.5) and {Gr = 0.9, Ga = 0.42) 
with e = 1.0 for the modeled micro-roughness that exhib- 
ited no hydrodynamic slip in Sec. IIII A\ and {Gr = 1.33, 
Ga = 0.57) with e = 1.3 for the surface that exhibited 
hydrodynamic slip in Sec. IIII Al Contact angle values are 
numerically measured via circle fitting within the hydro- 
dynamic region {(j)s < 0.01) with a procedure similar to 
that employed for the static contact angle. As showed in 
Fig. [9l advancing and receding angles are evaluated with 
the radii of the two circles that fit the rear and front 
interfaces connecting at height h. 

The applied pressure difference {Ap = 0.175-5.25 
xlO~^) causes a linear increase in the mean droplet ve- 
locity U = 0-0.01; the advancing and receding contact 
lines move with a velocity approximately equal to the 
computed mean speed U. Thus we determine a capillary 
number Ca = UfiL/j which we employ in Fig. [10] to re- 
port the measured dynamic contact angles. The key ob- 
servations from Fig. [To] are the following. (1) The model 
produces a small but finite contact angle hysteresis in 
static conditions {Ca 0) for a macroscopically smooth 
surface. We remark that we only have considered rough- 
ness at a microscale comparable the range of action of 
surface forces, which we modeled with a wall probability 
(2) The static hysteresis of the contact angles is pro- 
duced by the functional shape of the disjoining pressure, 
as predicted in |9|, and shows no significant dependence 
on the modeled microscale roughness. (3) Advancing and 
receding contact angles vary at different rates; while the 
advancing contact angle 6a increases with Ca, the reced- 
ing contact angle Or is weakly affected. This last obser- 




^^o-^r "^Gr = 1.33 Ga = 0.5 (e = 1.0) 
-Q-GR = lMGA=0S7{e = 1.3) 
■v Gr = 0.9 Ga = 0.42 (e = 1.0) 



0.005 0.01 0.015 0.02 0.025 0.03 
capillary number Ca 



FIG. 10. Advancing and receding contact angles versus cap- 
illary number Ga = U/jl/j (for the definition of the appro- 
priate U see text). All sets of interaction parameters result 
in the same value for the equilibrium angle Oy — 58° but 
different contact angle dynamics. The functional form of the 
disjoining pressure produces a finite static hysteresis of the 
contact angles. 



vation indicates the influence of the flow structure in the 
vicinity of a moving contact line and is further discussed 
in the next Section. 



D. Flow structure near a contact line 

The flow field in the vicinity of the advancing and re- 
ceding contact hues is reported in Figs. [TT1IT21 There 
are several key features produced by the present model 
whose existence and crucial exhibited effects have been 
discussed in previous works 0, H H. At the ap- 
parent three-phase contact line the flow becomes more 
intense, while the hydrodynamic slip significantly in- 
creases. Within the droplet bulk, away from the contact 
hues, little or no slip is observed at the hydrodynamic 
boundary (where (j)s{yw) = 0.01). The flow kinemat- 
ics approaching the advancing contact line resembles the 
corner flow configuration assumed by classical Voinov- 
Cox models [3] . 

The flow at the receding contact line, however, de- 
velops a trailing vortex [see Fig. [T2]f b-c)] as the droplet 
moves forward. Similar trailing vortices have been exper- 
imentally reported using PIV visualization in [26|. The 
reported flow features elucidate the widely different dy- 
namics of the advancing and receding contact angles re- 
ported in Fig. [To] for the modeled volatile fluid on macro- 
scopically smooth surfaces. The substantial difference in 
the flow kinematics at the receding and advancing con- 
tact lines can be thought to result from the combined 
effects of a large gradient of the disjoining pressure and 
hydrodynamics in the droplet bulk. 

We close this section with a few additional remarks. 
It is widely accepted that some form of hydrodynamic 
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slip must occur in order to allow the motion of a contact 
line ^-|6j. Previous works [H, [Bl, Q pointed out, however, 
that simply replacing no-slip by a slip boundary condi- 
tion does not preserve the experimentally observed kine- 
matics of the flow. Physical and conceptual issues remain 
in the modeling of hydrodynamic boundary conditions for 
moving contact lines [see [2] for a discussion]. Further- 
more, there is some experimental [24] and theoretical [2^ 
evidence that dynamic contact angles cannot be solely a 
function of the contact line speed and the physical prop- 
erties of the system. A mesoscopic model does not re- 
quire prescribing hydrodynamic boundary conditions nor 
the value of dynamic contact angles and, thus constitutes 
an interesting alternative to study dynamic wetting phe- 
nomena. Furthermore, conditions within the solid phase 
are dynamically computed with the proposed mesoscopic 
model and one does not need to prescribe any variables 
at the fluid-solid interface (e.g. pseudo-potential value, 
gradient of an order parameter), as it is customary for 
other LB models 0,0. 



IV. SUMMARY 

We have presented a mesoscopic model, based on a 
statistical description of the microscopic physics, that is 
capable of simulating nontrivial macroscopic phenomena 
observed during dynamic wetting of solids when hydro- 
dynamic effects play a critical role. Similarly to pre- 
vious models based on the lattice Boltzmann method 
[l0|, [ill , long-range microscopic interactions are modeled 
via mean-field potentials (i.e. pseudo-potentials) that are 
explicit functions of macroscopic variables. The most 
distinct feature of this model is that the solid phase is 
treated on a similar footing as the fluid phase. Short- 
range fluid-solid interactions are modeled by a solid prob- 
ability function (ps-, defining the probability of fluid-solid 
collisions, while long-range fluid-solid interactions are de- 
termined by pseudo-potentials designed to produce a dis- 
joining pressure with the proper functional form. Hence, 
the modeled surface forces have essential features re- 
quired to mimic physical effects observed during both 
full and partial wetting of solids. Regarding these funda- 
mental effects we remark the following. (1) The effective 
slip velocity, determined by momentum fluxes across the 
surface film, can widely vary on surfaces that exhibit the 
same equilibrium contact angle (i.e. surfaces with similar 
chemical properties). (2) The magnitude of the hydro- 
dynamic slip can significantly vary across macroscopic 
scales; the slip velocity increases near a moving contact 
line, while no slip is approximately recovered far from it. 
(3) Precursor films develop on hydrophilic surfaces and 
their evolution is dictated by the EOS of the fluid and 
surface forces; for the simulated system a small fraction 
of vapor is adsorbed on the solid substrate underneath 
the film. (4) A finite hysteresis of the static contact an- 
gle is observed for partially wettable and macroscopically 
smooth surfaces; the magnitude of this hysteresis is de- 




Os=0.5 



(b) 

FIG. 11. Fluid momentum magnitude IM^^-*! (top panels in 
the subfigures) and velocity streamlines (bottom panel in the 
subfigures) for droplet motion under pressure-driven flow, (a) 
Static contact angle By — 58° {Gr — 1.33, Ga — 0.5, e = 1.0). 
(b) Static contact angle By = 113° {Or = 1.33, Ga = 0.2, 
e= 1.0) 



termined by the modeled shape of the disjoining pres- 
sure. (5) The flow kinematics near advancing and re- 
ceding contact lines can exhibit substantial differences; 
rolling motion is favored near the receding contact line, 
while corner flow easily develops near the advancing con- 
tact Hue. The ability of the model we presented to re- 
produce these physical phenomena is crucial for studying 
diverse problems involving partial wetting and hydrody- 
namic interactions. 
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(d) 

FIG. 12. Flow patterns in the vicinity of the advancing and 
receding contact hnes observed on a hydrophihc surface and 
on a hydrophobic surface, (a) Advancing contact hne, static 
contact angle By = 58° {Gr = 1.33, Ga = 0.5, e = 1.0). (b) 
Receding contact line, static contact angle By = 58° {Gr = 
1.33, Ga = 0.5, e = 1.0). (c) Advancing contact line, static 
contact angle By = 113° {Gr = 1.33, Ga = 0.2, e = 1.0). 
(d) Receding contact line, static contact angle By — 113° 
(Gi? = 1.33, 6^^ = 0.2, e=1.0). 



V. DISCUSSION AND CONCLUSIONS 

There are several limitations inherent in the LB imple- 
mentation that originate with the discretization of phase 
space. As a consequence, the method is stable and ro- 
bust within a range of moderate density and viscosity 
ratios between the vapor and liquid phases, while it is 
difficult to realize large surface tensions (i.e. above one 
in lattice units). The study of many systems of inter- 
est can be accomplished by carefully matching the most 
relevant dimensionless groups in the problem of interest 
(e.g. Capillary number, Bond number, Ohnesorge num- 
ber). However, modeling of systems with large density 
ratios (e.g. water-air) in general flow conditions remains 
a challenging task and requires further developments of 
LB and other mesoscopic methods. An interesting possi- 
bility is to implement a mesoscopic approach at the same 
level of modeling by employing a particle-based method; 
this could potentially extend the range of physical pa- 
rameters where the numerical method is stable at the 
expense of certain computational advantages in the lat- 
tice discretization. 

We should remark that the proposed approach aims 
to reproduce the proper macroscopic physics in the fluid 
bulk. The purpose of our mesoscopic description is to dy- 
namically couple surface forces and hydrodynamics, cir- 
cumventing the use of complex boundary conditions for 
macroscopic variables. The inter facial film thickness in 
the proposed model is arbitrarily prescribed by placing 
a minimal number of simulation nodes (i.e. one lattice 
cell) required to resolve the interfacial region [see Ap- 
pendix]. The most attractive property of our approach 
lies in the fact that the chemical physics and hydrody- 
namics, fully coupled in the mesoscopic model, arise as 
different aspects of the simulated dynamics of a single- 
particle distribution. This level of description has the 
potential of guiding the study of complex phenomena in 
contact line motion. Such phenomena could include the 
finite time for formation and relaxation of the solid-fluid 
interface, and fluid compressibility effects. The applica- 
tion of this model to study deformable solids is of partic- 
ular interest; there, (j)s becomes time dependent. Other 
interesting applications include the dynamic wetting of 
micro-structured surfaces, and the extension of the pro- 
posed model to multi-component systems (e.g. colloidal 
suspensions, complex fluids). 
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Appendix A: The wall functions 

All fluid-solid interactions are determined by the 
spatially- varying functions (05, {pR^ and {^a) shown in 
Fig. [21 which must be defined before the dynamic sim- 
ulation starts. In principle, the three functions can be 
any integrable function with a finite value in the solid 
bulk and null value within the fiuid bulk. The spatial 
variation of these functions determines the magnitude of 
the surface forces (i.e. disjoining pressure) and momen- 
tum fiuxes in the dynamic simulation. The wall func- 
tion 05 (x) introduced in Eq. [12] controls the momen- 
tum exchange at the fiuid-solid interface and thus the 
effective slip velocity dynamically attained above the in- 
terfacial film. A finite value of the local wall probabil- 
ity (0 < 05 < 1) can be interpreted as the presence of 
microscopic-scale imperfections (or sub-nanometer scale 
roughness) on the solid surface; the product p05(x) rep- 
resents the local probability of a collisional event (short- 
range repulsion) between solid and fiuid molecules. The 
repulsive component, t/J^^, and attractive component 



of the Fluid-Solid pseudo-potential, ^I^fs-, give rise to a 
disjoining pressure with both repulsive and attractive 
parts; this property is crucial to model partially wettable 
surfaces [9]. 

The conditions described above could be easily satis- 
fied by simple analytical expressions in the case of fiat 
surfaces, as we studied in this work. However, it is con- 
venient to have a general numerical procedure to deal 
with complex surface geometries. For that purpose we 
employ a recursive Gaussian filter 

A^. Q 

Gf'\^) = ^ ^t..4--^)(x + r,). (Al) 

The initial wall probability = i^(x) is a Heaviside 
step function that takes the unit value, i^(x g 1^) = 1, 
in the region containing the solid nodes and takes the 
zero value outside this region, i^(x ^ 1^) = 0. In our sim- 
ulations the outer boundary of the domain dVt belongs to 
the solid bulk, in this region the value of the wall prob- 
ability must always be unit G^g^'^dn = 1- Accordingly, 

the smoothed wall function G^^^""^ is reset to unity on the 
outer boundary nodes dft after each iteration in Eq. IA1[ 
All numerical results in this work employ tpR = , 
= Gf^ - Gf, <Ps = Gf^ [see Fig.|2]. 
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